####################################################################
# File-Name: ItPSR simulations for graph                           #
# Date: November 2015                                              #
# Author: Ruth Dassonneville                                       #                      
####################################################################

## set working directory [change working directory if necessary]
setwd("C:/Documents and Settings/ruth dassonneville/Bureaublad/")

## load packages
library(foreign)
library(plm)
library(lme4)

## read data [adjust path if necessary]
data <- read.dta("/Users/ruthdassonneville/Desktop/ItPSR_Punishing Local Incumbents_replication.dta")


#############################################################################################

## Logit model explaining voting for incumbent

M1 <- glmer(IncumbencyCoalition ~  LaggedCoalition + Age+ Female + French +
  HighEdu+ Years_municipality+ Left_right + Terms + ENP + Unempl_change_10_11
           + (1|Municipality_code), family="binomial", data=data) 
summary(M1)

##get coefficients and varcov matrix
est<-fixef(M1)
V<-vcov(M1)
est
V

## set up sampling distribution
library(MASS)
N<-10000
S<-mvrnorm(10000, est,V)
head(S)

##choose sensible covariates
means<-apply(model.matrix(M1), 2, mean)
means

## predictions Unempl -2 +10 - [adjust several times for each point prediction]

X2<-means
X2
X2[11]<-10
head(X2)

## get linear predictor and form expected value
mu<-S%*%X2
mean(mu)
hist(mu)

lower<-apply(mu,2,quantile, 0.2)
lower
upper<-apply(mu,2,quantile, 0.8)
upper

#in response function to get predicted probability
p<- 1/(1+exp(-mu))
mean(p)

plower<-1/(1+exp(-lower))
pupper<-1/(1+exp(-upper))
plower
pupper

#############################################################################################

